Glass transition and effective potential in the hypernetted chain approximation 
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l~^ ■ We study the glassy transition for simple liquids in the hypernetted chain (HNC) approximation 

O^ ' by means of an effective potential recently introduced. Integrating the HNC equations for hard 

spheres, we find a transition scenario analogous to that of the long range disordered systems with 
_i "one step replica symmetry breaking" . Our result agree qualitatively with Monte Carlo simulations 

"T^ I of three dimensional hard spheres. 

D 



"^ . The HNC approximation is one of the most widely used approach to describe the density-density correlation function 

g{x) for hquids at equilibrium ||l|. It consists in a self-consistent integral equation that can be derived by a partial 
C^ , resummation of the Meyer expansion, and corresponds to the variational equation for a suitable free-energy functional 
P,||. The simple HNC approach does not by itself allow to detect freezing M. The simple inspection of the pair 
correlation function certainly does not allow to do so, being qualitatively similar in the liquid and in the glass. 
■ O ■ Freezing, although present, can be hidden if one concentrates on simple equilibrium quantities B. 

In has been recently stressed in 1^ that the freezing transition can be detected combining the HNC approximation 
with the replica method by studying the correlation functions among different replicas of the same system in presence 
of a potential which couples them. At low temperature (or at high density) one finds a self-consistent solution where 
different replicas remains correlated also in the limit of zero coupling. This phenomenon correspond to freezing and 
it goes under the technical name of replica symmetry breaking. 

In this letter we pursue this idea of studying the glass transition in the HCN approximation. We are not concerned 
Qv' ' about the behavior in the glassy phase. Our aim it to use an effective potential recently introduced by two of us [^,0, 
f^ to study the glass transition of HNC hard spheres in 3D. We compare the results with Monte Carlo simulations of real 
P\| I hard spheres. The conceptual advantage of this approach is that all the subtle points of the usual approach related 
^~H to replica symmetry breaking are not needed in order to expose the transition. 

r^ . The effective potential is constructed as follows, for a system described by the coordinates of all the particles 

^^ ] X = (xi, ...,a;jv) and with potential energy H{x) = X)i<i 4'{xi ~ ^j)- 

jrt ■ Let us consider a reference configuration y chosen with probability exp(— /3'iJ(j/))/Z(/3'), where /?' = 1/T' is 

Cj [ some arbitrary inverse temperature. Let us define a distance among configuration as d{x,y) = 1 — q{x,y), with 

I ■ the "overlap" q{x,y) defined as q{x,y) = jfj^i'i w{\xi — yj\). w is an attractive potential that in this letter we 

^ , choose as w{r) — 9{rQ — r). To very different configurations it corresponds large distance and small overlap, to similar 

configurations small distance and large overlap. We define a constrained Boltzmann-Gibbs measure at temperature 

T as 
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KAV) = yi^ . e-^^(-) 5{q{x, y) q) (1) 

Z[P,q,y) 



j^ • where Z{f3,q,y) is the integral over x of the numerator. This conditional measure allows to probe regions of the 
configuration space having vanishingly small probability, and as we will see, it will help us to reveal the glassy 
structure hidden in the simple equilibrium approach. 

Introducing a Lagrange multiplier conjugated to q to enforce the delta function and integrating over it by saddle 
point, one see that the free energy associated to (|^), V{q) = —T\ogZ{P,q,y), can be computed as the Legendre 
transform of F(e) = — T log Z{P,e,y) with Z{(3,e,y) = J dx e~^^-^^^^~'^'^^^'^''\ If the coupling e is positive there is 
an attraction to the reference configuration y. Of special interest will be the cases e — > 0"*", while q will go to a non 
trivial value. The free-energy F(e) and the potential V{q) should be self-averaging with respect to the distribution 
of y, and be therefore just functions of their argument and the temperatures /3 and /?'. From now on, in this letter 



we will limit ourselves to the case j3 = (3' which will be enough to detect freezing in the system. It is conceptually 
important however to consider the more general case if one would like to describe a system which, after crossing the 
freezing temperature, remains confined in the vicinity of the configuration where it was last able to thermalize. 
In order to compute F(e) in any physical system we need to average Z[(3, e, y) over the distribution of y. This can 

be done in a convenient way using the replica method, where one writes \ogZ = lim^^o , and compute the limit 

from an analytic continuation from integer r. In principle the replica method can be avoided but it is quite useful to 
make all the computation quite straightforward. Explicitly: 



Z'' = / da::oda:i...da;r e 



-'^ EI=o "i^<^)+p^ EI=i li^o.x^ 



(2) 



we have written xq = y. The problem is reduced to that of an equilibrium mixture of r + 1 species (with r — > 0), 
and is formally similar to the one developed by Given, StcU and collaborators to study liquids in random matrices 
[p|. The use of the formalism is however different. In S the replica method was used to deal with the quenched 
disorder represented by the medium, while for us the potential is a tool to probe regions of configuration space of 
small Boltzmann probability and_ we do not have quenched disorder. The HNC equation can be derived from the 
following free-energy functional 



WF{e) 



A'^x 



E 

a, 6=0 



P^gab{x) [\0ggab{x) - 1 + I3(j){x)5ab\ 



2/3e^p^goa(a;)w(a;) 



Tr L(/9/i) 



(3) 



with L(?i) = u — V? /2 — log(l + u), hab = fjab ^ Ij the trace of L is intended both on replica indices and in the operator 
sense in space. The equations, and the relative value of the free energy are obtained extremizing (0) over all the 
replica correlation functions Safe's, and extracting the terms of order r. 
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FIG. 1. The behavior of g as a function of e for HNC hard spheres for p — 1.14, 1.17, 1.19, 1.20. For high enough density q 
is a multivalued function of e. We have shown only a portion of the curve in the region where it is multivalued. For graphical 
transparency in this and the next figure we have joined with a line the branches corresponding to the same density. 



In order to continue analytically F we use the ansatz gab — goo for a — b — 0, gab = gio for a = or 6 = and 
a ^ h and gab = 5^6 fo'" both a and b different from 0. In this paper we will only consider the replica symmetric 
ansatz g*^ = gii5ab + 512(1 — Sab)- We warn the reader that this ansatz should give the correct value of the potential 
V{q) for high and low values of q in the liquid phase, replica symmetry breaking effects are to be expected in an 
intermediate regime H even in the liquid phase. The physical meaning of the various elements of gab within this 
ansatz is immediate. The element goo represents the pair correlation function of the free system; as such the equation 



determining it decouples from the other components in the Hmit r — > and coincides with the usual HNC equation 
for a single component system. In turn, gn represents the pair correlation function of the coupled system. giQ is 
the pair correlation among the quenched configuration and the annealed one, while 1712 represents the correlation 
between two systems coupled with the same quenched system. This last is the analogous of the Edwards- Anderson 
order parameter in disordered systems pOf , and represents the long time limit of the time dependent autocorrelation 
function at equilibrium te[. 
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FIG. 2. The effective potential for HNC hard spheres. From to to bottom p = 1.0, 1.14, 1.17, 1.19, 1.20. For low density, 
high up in the liquid phase the potential is convex. In the glass phase two minima are present. 



In order to study g as a function of e we have solved the HNC equations in 3D with the hard sphere potential 
4>{r) = 00 for r < 1, <j){r) = for r > 1 for various values of the density (using a space resolution equal to .01 and a 
large distance cutoff equal to 10). The density p is the control parameter of the freezing transition in this problem 
where there is no temperature. We have rescaled /3e — > e, and chosen rg = 0.3 in the definition of the overlap. We 
have reconstructed the curves q{e) and V{q) following the solution of the HNC equations starting from high and low 
values of e and respectively decreasing or increasing it slowly. We see in figure |l] and g that for low enough density 
g is a single valued function and the potential is convex, with the minimum corresponding to 501(2^) = 1 for all x. 
This is a fair description of the liquid phase. Above a critical density p* « 1.14 the potential loose convexity, and a 
coupling can induce a transition between a high and a low q phase H. For p = pc ^ 1.17 a second minimum at high 
q appears, deepening and deepening as the density is increased. The presence of this second minimum shows that 
above pc the system with is in the glassy phase. If, by means of a large e we prepare the system in the vicinity of y 
and we then let e -^ 0, the system remains confined. It should be noticed that while the shape of V{q) depends on 
the particular definition of the overlap, the properties of the minima of the potential do not, as they correspond to 
vanishing coupling. 

As it has been discussed in [HQ], the two minima structure is associated to a Gibbs-Di Marzio glass transition 
scenario. At pc the ergodicity is broken and an extensively large number of metastable states Af = e^^ contributes 
to the partition function. The relative height of the two minima is exactly equal to the "configurational entropy" E, 
and as we can see in fig. pi as p is increased E decreases, until it vanishes at p = ps « 1.203. (The values of pc and 
Pc are compatible with those found in 0], indeed the potential method reproduces the results of replica symmetry 
braking approach for the static and dynamic critical densities). The shape of the potential is the characteristic one of 
a system undergoing a first order phase transition. We can use Maxwell construction to locate the transition line in 
the plane e — p, which is shown in fig. |4| The computation as it stands is not consistent for p > ps'. it gives a negative 
configurational entropy in that region. To describe consistently the behavior there the replica symmetry breaking 
formalism of |3|| is needed. 
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FIG. 3. The configurational entropy E as a function of p. 

Although we did not try to compare quantitatively the values of the freezing density we get with the one previously 
found in numerical simulations [|ll| , we have performed our own Monte Carlo simulations to test in a qualitative way 
the prediction of a first order transition in presence of a coupling with a fixed configuration. 
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FIG. 4. Phase diagram in the plane e — p. A first order transition line terminating in a critical point separates a low q from 
a high q phase. 



To generate configurations at fixed density we start with N particles of zero radius in a three dimensional box with 
periodic boundary conditions, and we let the radii grow until two particle do get in contact. At that point we make a 
Monte Carlo sweep and iterate the procedure until the desired density is reached. The volume and the radius (r) are 



at the point rescaled in order to have r = 1. We thermaUze then the system for 4000 Monte Carlo sweeps and use the 
configuration y reached as "external field" for our coupled replicas experiments. The relatively short thernialization 
is chosen in order to avoid crystallization. Having generated the configuration y we let evolve a coupled system x. For 
various densities, we start the evolution from the configuration y with an high value of e and decrease the value of e in 
units of (5e, making 2*^ Monte Carlo iterations for each value of e. In figure p^wc plot g as a function of e for different 
values of k. We see that, as it should be expected for a system undergoing a first order phase transition, the curves 
are smooth for low k and tend to develop a discontinuity for large k. We have presented here results for p = 1.04. 
Other simulations (which we do not display here) show lower density the discontinuity occurs at higher e, while it is 
pushed toward smaller e for higher density. From the quantitative point of view there is about a 20% agreement on 
the value of the density at which a transition is present in the e-q plane; however the qualitative prediction of a first 
order transition at vales of e of order 1 is clearly satisfied. 




FIG. 5. The behavior of g as a function of e for a system of 258 particles and p — 1.04. The different curves correspond 
different thernialization times 2* for each value of e. ^From top to bottom k — 17, 19, 21, 23. For larger thermalization times 
the system seem to develop a first order jump in q. 



A different numerical experiment is presented in figure H. Here we let the system evolve at fixed e starting at time 
zero from x = y and we plot the overlap as a function of time. Again we observe a behavior compatible with a 
discontinuity of g as a function of e. 
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FIG. 6. The overlap as a function of time in a logarithmic scale (the horizontal axis is ln2(i) starting from y at time and 
evolving for fixed e. In this figure p = 1.04 and from top to bottom e — 0.8, 0.7, 0.6, 0.5, the number of particles is 1024. 

We see that the HNC approach predicts a glass transition scenario very close to the one found in systems with 
"one step replica symmetry breaking" with a non convex |g] effective potential and Gibbs-Di Marzio entropy crisis 
| p2[ . The HNC is in this respect a genuine mean-field theory, giving infinite life meta-stable states. In real systems 
mctastable states have finite life time and the potential has to be a convex function of q for any densities. As it has 
been discussed many times in the application of the theory to real systems [Q the picture should be corrected to take 
that into account. The density pc, representing the point where the relaxation time diverges in mean-field, becomes 
a crossover value where the dynamics is dominated by barrier jumping processes [y. Elucidation of the dynamical 
processes responsible for restoration of ergodicity beyond mean-field is one of the currently open issue in glass physics. 
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